Hyper-sparsity of the density matrix in a wavelet representation 
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O(N) methods are based on the decay properties of the density matrix in real space, an effect 
sometimes refered to as near-sightedness. We show, that in addition to this near-sightedness in 
real space there is also a near-sightedness in Fourier space. Using a basis set with good localiza- 
tion properties in both real and Fourier space such as wavelets, one can exploit both localization 
properties to obtain a density matrix which exhibits additional sparseness properties compared to 
the scenario where one has a basis set with real space localization only. We will call this additional 
sparsity hyper-sparsity. Taking advantage of this hyper-sparsity, it is possible to represent very large 
quantum mechanical systems in a highly compact way. This can be done both for insulating and 
metallic systems and for arbitrarily accurate basis sets. We expect that hyper-sparsity will pave the 
way for O(N) calculations of large systems requiring many basis functions per atom, such as Density 
Functional calculations. 

PACS numbers: 71.15-m 



Methods for the calculation of the electronic structure 
that exhibit linear scaling with respect to the size of the 
system, so-called O(N) methods Q, are an important 
topic in physics and chemistry. Only with this kind of 
algorithms it is possible to calculate the properties of 
very large systems that are relevant in many practical 
applications. O(N) methods have become practically a 
standard for tight binding calculations where the ratio 
of the number of basis functions to the number of elec- 
trons is of the order of 2. In highly accurate Density 
Functional type calculations, where the number of basis 
functions is very large compared to the number of elec- 
trons, these methods have on the other hand not been 
widely used up to now. This comes form the fact that 
O(N) methods are only efficient if the density matrix is 
very sparse. Large conventional basis sets lead however 
to density matrices that are not sufficiently sparse. The 
sparseness of the density matrix in connection with con- 
ventional localized basis sets is due to the real space decay 
properties of the density matrix, the so called real space 
near-sightedness Q of the system. We will show how an 
additional near-sightedness in Fourier space can be used 
to make the density matrix significantly more sparse. 

For reasons of simplicity and in order to be able to do 
certain numerical calculations without any truncation we 
will concentrate in the following on the one dimensional 
case even though all the principles are applicable to the 
three dimensional case. For our investigations we use a 
test system that is characterized by a simple harmonic 
potential —2 sin(27rr) and we applied periodic boundary 
conditions. The valence band extends from roughly to 
4 atomic units and it is followed by a rather large gap of 
roughly 2 atomic units. By occupying each primitive cell 
that has a length of one atomic unit with one electron 



pair one thus obtains an insulator, while one obtains a 
metal by assigning for instance 2 electron pairs to each 
cell. 

The density matrix F in a real space representation 
can be modelled by 

sinf— (r — r')) 
F(r,r') = C exp( K r) — ^ — (1) 

r — r 

This is in principle only the asymptotic form, but is also 
a good approximation for small values of r and r'. The 
decay constant k grows with increasing gap and is zero 
in a metallic system (3) . The oscillation length a is of the 
order of the average distance among the valence electrons 
in a metal and of the order of the interatomic spacing in 
an insulator. C is a normalization constant. Evidently 
F(r, r') has both a typical length scale of — and a domi- 
nating Fourier component of — . The localization in real 
space will be particularly good for strongly insulating sys- 
tems whereas the localization around the typical Fourier 
component becomes better for decreasing gap sizes. It is 
well know that in most realistic materials F(r, r') decays 
over several oscillation lengths indicating a combined lo- 
calization in real and Fourier space. 

Let us now briefly review some essential facts from 
wavelet theory Q . For simplicity we will work with the 
standard Daubechies wavelet family There are two 
sets of fundamental functions, the scaling functions <p and 
the wavelets ip. The scaling functions are essentially ordi- 
nary localized functions comparable for instance to Gaus- 
sians in many respects. Their advantage over Gaussians 
is however that they are orthogonal. To form a basis set 
one has to take all the integer translations i of the scaling 
functions at a certain resolution level k 

$(x) = cfr(2 k x-i) (2) 
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Applying the unitary fast wavelet transform, one can 
transform the set of the scaling functions <fi into the 
equivalent set of wavelets tp. The wavelet basis one ob- 
tains in this way consists now of translated wavelets 
at different resolution levels. What is important in 
this context is that the scaling functions do not ex- 
hibit any localization in Fourier space whereas the 
wavelets do exhibit this feature. The Fourier spec- 
trum for the Daubechies wavelet of degree 12 (Fig- 
ure 0) that was used in this work is shown in Figure [|. 
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FIG. 1. The Daubechies scaling function (solid line) and 
wavelet (dashed line) used in this work. 




FIG. 2. The Fourier power spectrum P(k) of the 
Daubechies wavelet at three different resolution levels denoted 
by skinny, normal and fat. As one sees the peaks at the domi- 
nating Fourier components are reasonably well separated from 
each other. 

Let us now go over to O(N) context. The most general 
O(N) methods are the ones that calculate the full den- 
sity matrix. In these methods one calculates the matrix 
elements of F with respect to a basis set \i- 



F, 



Xi(r)F(r, r')xj(r')drdr' 



(3) 



Let us first discuss the case where the x's are ordinary 
localized functions. In order to be able to do a fair eval- 
uation of the advantages obtained by additional Fourier 
localization we will chose as our set of localized functions 
the set of scaling functions that is completely equivalent 
to the set of wavelets used in the comparison. It is obvi- 
ous that the width of the scaling functions has to be less 



than the wavelength of the oscillatory part of F in order 
to be an accurate basis. This means that 



p. . — 



$(r)F(r,r r )tf(r')drdr' « F(Ri,Rj) (4) 



where Ri is the position of the z-th scaling function. The 
function -F(0, r) for our test system is shown in Figure [|. 
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FIG. 3. The function F(0,r). Along the left y-axis going 
with the solid line it is plotted on a normal scale and along 
the right y-axis going with the dashed line on a logarithmic 
scale. 

Inspite of the exponential decay, one has to allow for 
a fairly large sub-volume of the whole 64 atom sys- 
tem before one can truncate F(0,r) without a large 
truncation error. This means, that most of the ele- 
ments of the density matrix have to be calculated and 
that therefore any traditional O(N) scheme will not 
be very efficient. The structure of the density ma- 
trix in this case is shown in Figure ^. Four scal- 
ing functions per atom were used, which gives a to- 
tal energy that deviates by 4.e-3 a.u. per atom from 
the exact result in the limit of an infinite basis set. 
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FIG. 4. Structure of the density matrix in a scaling func- 
tion representation. Elements larger than l.e-3 are denoted 
by green areas elements larger than l.e-2 by blue areas. 

Let us now look at the same matrix in a wavelet basis, 
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where we have to calculate matrix elements of the type 

ipi(r)F(r,r')ipf (r')drdr' (5) 

Note that in contrast to Equation ^ there are now matrix 
elements between wavelets at different resolution levels k 
and fc', that have different dominating Fourier compo- 
nents. In accordance with the usual convention we order 
the wavelets at the different resolution levels such that 
the skinny wavelets follow the fat ones. So in our exam- 
ple, the indices 129 to 256 refer to the skinniest wavelets, 
the indices 65 to 128 to the not quite so skinny wavelets 
etc. Let us first look at the upper right 128 times 128 
block in Figure ||. This block represents only matrix el- 
ements among the most skinny wavelets. From Figure ^ 
we would expect to find elements larger than l.e-2 and 
that the bandwidth is close to half of the size of the 
block. In reality all elements are smaller than l.e-2, and 
the bandwidth is much less than the one in Figure ^. 
This is the effect of the above postulated hyper-sparsity. 
If both basis functions in Equation ^ are skinny wavelets, 
the function F(r, r') can not couple them since its dom- 
inating wave length is much larger than the dominat- 
ing wave length of the two wavelets. As a matter of 
fact any matrix element will be strongly damped un- 
less all the three terms in Equation |j| have compara- 
ble wave lengths. Therefore also all the other fingers 
in Figure ^ that represent coupling between wavelets of 
different stature (i.e. different k and k') have a small 
bandwidth. Because their spatial frequencies do not 
match, the decay for wavelets positioned a certain dis- 
tance apart is much faster than one would expect from 
the decay behavior of the amplitude of F(r,r') alone. 
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FIG. 5. Structure of the density matrix in a wavelet rep- 
resentation. The convention for the colors is the same as in 
Figure [|. The decay is so fast in this case that the green area 
around the blue area is so thin to be hardly visible in some 
regions. 

To quantitatively assess the advantages of the addi- 



tional Fourier localization, compared with spatial lo- 
calization only let us plot the number of coefficients 
that are necessary to represent the density matrix with 
a certain error. We define the error as the 2-norm 
of the difference between the exact and the truncated 
density matrix. As one can see from Figure ^ the 
effect is dramatic. The accuracy obtained with the 
same number of coefficients with the wavelet represen- 
tation is several orders of magnitude smaller than the 
one obtained with the scaling function representation. 
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FIG. 6. The error in the density matrix versus the size 
of the data set necessary for its representation. Solid lines 
correspond to a scaling function representation, dashed lines 
to a wavelet representation. The two red curves are for a 
metallic system, the two green curves for an insulator and 
the two blue curves for an insulator where the density matrix 
was constructed indirectly via the Wannier functions. 

In the case of an insulator at zero electronic tempera- 
ture, the density matrix can compactly be represented in 
terms of Wannier functions. The Wannier functions are 
characterized by a spatial decay length and a dominant 
Fourier component as well. So the same principles apply 
and a wavelet representation is significantly more efficient 
than a scaling function representation (Figure [sj) . Since 
one needs only one Wannier function per electron pair 
instead of all the columns of the density matrix the pref- 
actor in the Wannier representation is however smaller. 

Let us next consider the metallic systems. For the 
case of jellium, plane waves would obviously give the 
most compact representation, leading to a strictly diag- 
onal density matrix. So the principle of spatial near- 
sightedness does not apply and the Fourier space lo- 
calization of the basis functions will even be more im- 
portant than the real space localization in the metallic 
case. As one sees from Figure [l] our wavelets have a 
more pronounced spatial localization than Fourier local- 
ization. Other families of wavelet which have a more 
pronounced Fourier space localization might in this con- 
text be more efficient than the Daubechies wavelet used 
in these tests. The structure of the matrix is shown 
in Figure |?] for the case of the wavelet representation 
of the test system containing 128 electron pairs and 
512 basis functions. In the case of a scaling func- 
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tion representation all elements are larger than l.e-2 
and the corresponding matrix is therefore not shown. 
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FIG. 7. Same as Figure || but for a metallic system. 

As expected, the quantitative evaluation in Figure ^ 
shows that the savings of a wavelet representations com- 
pared to a scaling function representation are even more 
important in the metallic case than in the case of an 
insulator. In both representations more coefficients are 
however needed than in the insulating case. This comes 
from the fact that the decay behavior of matrix elements 
with moderate spatial frequency disparity is still deter- 
mined by the slow spatial decay of F(r, r'). In addition to 
the finger structures there are therefore some full blocks 
in Figure |?]. 

In the context of large systems, the scaling of the re- 
quired number of basis functions with respect to system 
size is of course important. As is well known, even with 
basis sets which have real space localization only, this 
scaling is linear. Hence it is to be expected to be linear 
as well in the case of a combined spatial and Fourier lo- 
calization. This assumption was clearly confirmed by our 
numerical experiments where we studied model systems 
containing between 16 and 128 " atoms" . For metallic sys- 
tems ordinary localized basis functions give a quadratic 
scaling unless one goes over to unrealistically large sys- 
tems. Our numerical test for metallic systems containing 
between 16 and 128 "atoms" gave an exponent of 1.66. 
So the advantages of basis sets with frequency localiza- 
tion will grow with increasing system sizes. 

Like in the majority of wavelet based calculations H[|, 
the basic computational kernels of the calculations pre- 
sented here were implemented with respect to a scaling 
function basis. The final results were then transformed 
into a wavelet basis to examine the influence of trun- 
cation. In order to take full advantage of a truncated 



wavelet basis set, it would be necessary to do the whole 
calculation within the truncated wavelet basis. Matrix 
times vector multiplications are the basic step of several 
O(N) schemes and it has been demonstrated how they 
can be implemented with strictly linear scaling with re- 
spect to the number of non-zero coefficients . So the 
basic observation of this paper could directly be used to 
reduce the numerical effort in the same way as the num- 
ber of necessary basis functions. 

Based on the observation that the density matrix 
F(r, r') is characterized by both a typical length scale 
and a dominating Fourier component, we have demon- 
strated that quantum mechanical systems exhibit not 
only a near-sightedness in real space but also a near- 
sightedness in Fourier space. A basis set that has com- 
bined localization in both real and Fourier space can 
therefore significantly reduce the number of data that are 
needed to describe it. We expect that this phenomenon 
that we call hyper-sparsity will allow to overcome the 
present limitation of O(N) scheme to systems where the 
number of basis functions per atom is rather small. We 
thank J. Flutter, M. Parrincllo, Anna Putrino and O. 
Jepsen for interesting discussions. 
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